// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                            avtH5PartFileFormat.C                          //
// ************************************************************************* //

#include <visit-config.h>
#include <avtH5PartFileFormat.h>

// VTK
#include <vtkCellData.h>
#include <vtkCellTypes.h>
#include <vtkUnsignedIntArray.h>
#include <vtkFloatArray.h>
#include <vtkDoubleArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>

// VisIt
#include <avtDatabaseMetaData.h>

#include <DBOptionsAttributes.h>
#include <Expression.h>
#include <avtCallback.h>
#include <InvalidVariableException.h>
#include <InvalidFilesException.h>
#include <NonCompliantFileException.h>
#include <BadIndexException.h>
#include <DebugStream.h>
#include <TimingsManager.h>

#ifdef PARALLEL
  #include <avtParallel.h>
#endif

// ****************************************************************************
//  Method: avtH5PartFileFormat constructor
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue eb 9 13:44:50 PST 2010
//
//   Modifications:
//    Gunther H. Weber, Thu Apr  1 18:17:11 PDT 2010
//    Enable domain decomposition per default.
//
//    Gunther H. Weber, Wed Aug  4 15:34:04 PDT 2010
//    Read and store number of time steps in file for use by GetNTimesteps().
//    Use lazy evaluation for queries. Do not run query in
//    RegisterDataSelections(). Instead, use flag "queryResultsValid" to
//    indicate whether current queryResults reflect current query.
//
//    Jeremy Meredith, Wed Mar 30 12:55:29 EDT 2011
//    Close the file if it's no an H5Part file -- HDF5 can be sensitive about
//    not closing files.
//
//    Gunther H. Weber, Wed Aug 17 16:28:46 PDT 2011
//    Write to debug log file if reader was compiled with FastBit support.
//
//    Eddie Rusu, Wed Oct 14 15:52:13 PDT 2020
//    Removed FastBit and FastQuery
//
// ****************************************************************************

avtH5PartFileFormat::avtH5PartFileFormat(const char *filename,
        const DBOptionsAttributes *readOpts) : avtMTSDFileFormat(&filename, 1)
{
    // ibis::gVerbose = 1000;
    int t1 = visitTimer->StartTimer();

    // Defaults
    enableDomainDecomposition = true;
    variablePathPrefix = std::string("Step#");

    // Check options
    if (readOpts != NULL)
    {
        for (int i = 0; i < readOpts->GetNumberOfOptions(); ++i)
        {
            if (readOpts->GetName(i) == "Enable domain decomposition")
                enableDomainDecomposition =
                  readOpts->GetBool("Enable domain decomposition");
            else if (readOpts->GetName(i) == "Variable path prefix")
                variablePathPrefix =
                  readOpts->GetString("Variable path prefix");   
        }
    }

    // only open the file briefly to ensure that it is really an
    // H5Part file
    file = H5PartOpenFile(filename, H5PART_READ);

    if (!file)
        EXCEPTION1(InvalidFilesException, filename);

    // This function was once supported by H5Part and will be again in the new
    // release. Until it is widely available, comment it out.
    if (H5PartFileIsValid(file) != H5PART_SUCCESS)
    {
        debug1 << "avtH5PartFileFormat::avtH5PartFileFormat(): "
               << "H5PartFileIsValid check failed." << std::endl;

        H5PartCloseFile(file);
        EXCEPTION1(InvalidFilesException, filename);
    }

    // Check for the first time step - if it does not exit bailout as
    // the file is not a H5Part file.
    activeTimeStep = 0;
    if (!H5PartHasStep(file, activeTimeStep))
    {
        H5PartCloseFile(file);
        EXCEPTION1(InvalidFilesException, "Cannot find first time step 0.");
    }

    // Assume a valid H5Part file from here out.
    
    // Get number of time steps in file
    numTimestepsInFile = H5PartGetNumSteps(file);

    // The following should be equivalent to testing the presence of a
    // group named "Group #0".
    if (numTimestepsInFile <= 0)
    {
        debug1 << "avtH5PartFileFormat::avtH5PartFileFormat(): "
               << "File contains " << numTimestepsInFile
               << " <= 0 time steps." << std::endl;

        H5PartCloseFile(file);
        EXCEPTION1(InvalidFilesException, filename);
    }

    for( int i=0; i<numTimestepsInFile; ++i )
    {
      if (H5PartSetStep(file, i) != H5PART_SUCCESS)
      {
          H5PartCloseFile(file);
          EXCEPTION1(InvalidFilesException, "Cannot activate time step.");
      }
      
      cycles.push_back( i );

      h5part_int64_t numStepAttr = H5PartGetNumStepAttribs( file );

      for( h5part_int64_t i=0; i<numStepAttr; ++i )
      {
        char attrib_name[maxVarNameLen];
        h5part_int64_t len_of_attrib_name;
        h5part_int64_t attrib_type;
        h5part_int64_t attrib_nelem;

        H5PartGetStepAttribInfo( file, i, attrib_name, maxVarNameLen,
                                 &attrib_type, &attrib_nelem );
        
        if( strncmp( attrib_name, "Time", 4 ) == 0 ||
            strncmp( attrib_name, "time", 4 ) == 0 )
        {
          h5part_int64_t status = H5PART_SUCCESS;

          if (attrib_type == H5PART_INT64)
          {
            h5part_int64_t time;
            status = H5PartReadStepAttrib(file, attrib_name, &time);
            times.push_back( time );
          }
          else if (attrib_type == H5PART_INT32)
          {
            h5part_int32_t time;
            status = H5PartReadStepAttrib(file, attrib_name, &time);
            times.push_back( time );
          }
          else if (attrib_type == H5PART_FLOAT64)
          {
            h5part_float64_t time;
            status = H5PartReadStepAttrib(file, attrib_name, &time);
            times.push_back( time );
          }
          else if (attrib_type == H5PART_FLOAT32)
          {
            h5part_float32_t time;
            status = H5PartReadStepAttrib(file, attrib_name, &time);
            times.push_back( time );
          }
          else
          {
            EXCEPTION2(NonCompliantFileException, "H5Part Open",
                       "Unsupported value type for coordinates. "
                       "(Supported tyes are FLOAT64, FLOAT32 and INT64.)");
          }
          // FIXME: Possibly handle other data types for coordinates.

          if (status != H5PART_SUCCESS)
          {
            EXCEPTION2(NonCompliantFileException, "H5Part Open",
                       "Could not read time.");
          }
        }
      }
    }
    
    // Activate first time step
    activeTimeStep = 0;
    if (H5PartSetStep(file, activeTimeStep) != H5PART_SUCCESS)
    {
        H5PartCloseFile(file);
        EXCEPTION1(InvalidFilesException, "Cannot activate time step 0.");
    }
  
    // Iterate over particle var names
    h5part_int64_t nPointVars = H5PartGetNumDatasets(file);
    if (nPointVars < 0)
    {
        H5PartCloseFile(file);
        EXCEPTION2(NonCompliantFileException, "H5Part Constructor",
                "Could not obtain the number of particle variables.");
    }

    // Get information about particle variables
    for (int i=0; i<nPointVars; ++i)
    {
        // Get information. We are interested in name and type
        char varName[maxVarNameLen];
        h5part_int64_t type;
        h5part_int64_t nElem;

        if (H5PartGetDatasetInfo(file, i, varName, maxVarNameLen, &type, &nElem)
                != H5PART_SUCCESS)
        {
            H5PartCloseFile(file);
            EXCEPTION2(NonCompliantFileException, "H5Part Constructor",
                       "Could not obtain particle data set information.");
        }

        // Store information about type in map so that we can access
        // it in GetVar()
        if (particleVarNameToTypeMap.find(varName) != particleVarNameToTypeMap.end())
        {
            H5PartCloseFile(file);
            EXCEPTION2(NonCompliantFileException, "H5Part Constructor",
                       "Particle data contains two variables with the same name.");
        }
        else
        {
            particleVarNameToTypeMap[varName] = type;
        }
    }

    particleNSpatialDims = 2;
    unsigned char cartesian = 0;
    unsigned char cylindrical = 0;
    unsigned char spherical = 0;

    // If there is particle data in the file (check needed since a
    // file could contain only block data).
    if (particleVarNameToTypeMap.size())
    {
        // Determine the coordinate system used and whether the file
        // is 2D or 3D
        for (VarNameToInt64Map_t::const_iterator it = particleVarNameToTypeMap.begin();
             it != particleVarNameToTypeMap.end(); ++it)
        {
            std::string currVarName = it->first; 

            std::transform(currVarName.begin(), currVarName.end(),
                           currVarName.begin(), ::tolower);

            if (currVarName == "x" || currVarName == "y" || currVarName == "z" )
                ++cartesian;

            if (currVarName == "r" || currVarName == "phi" || currVarName == "z")
                ++cylindrical;

            if (currVarName == "r" || currVarName == "phi" || currVarName == "theta" )
                ++spherical;

            if (currVarName == "z" || currVarName == "theta" )
                particleNSpatialDims = 3;

            
        }

        // Spherical coordinates always requires three values. 
        if( spherical == 3 && (cylindrical >= 3 || cartesian >= 3) )
        {
          std::stringstream buf;

          buf << "Ambiguous coordinate system: both "
              << "spherical and cylindrical or cartesian coordinate systems "
              << " were found. Can not disambiguate. ";

          debug1 << "avtH5PartFileFormat::avtH5PartFileFormat(): "
                 << buf.str() << std::endl;

          EXCEPTION2(NonCompliantFileException,
                     "H5Part Constructor", buf.str().c_str() );

          // callback is not working so use the exception above
          avtCallback::IssueWarning( buf.str().c_str() );
        }

        // Cylindrical and cartesian coordinates have two or three values.
        else if( cylindrical >= 2 && cartesian >= 2 )
        {
          std::stringstream buf;

          buf << "Ambiguous coordinate system: both "
              << "cylindrical and cartesian coordinate systems were found. "
              << "Can not disambiguate. ";

          debug1 << "avtH5PartFileFormat::avtH5PartFileFormat(): "
                 << buf.str() << std::endl;

          EXCEPTION2(NonCompliantFileException,
                     "H5Part Constructor", buf.str().c_str());

          // callback is not working so use the exception above
          avtCallback::IssueWarning( buf.str().c_str() );
        }

        if( spherical == 3)           { coordType = sphericalCoordSystem; }
        else if( cylindrical >= 2)    { coordType = cylindricalCoordSystem; }
        else /* if( cartesian >= 2)*/ { coordType = cartesianCoordSystem; }

        debug1 << "File contains a " << particleNSpatialDims
               << "D mesh using a "
               << (coordType == cartesianCoordSystem ? "Cartesian" :
                   (coordType == cylindricalCoordSystem ? "Cylindrical" :
                    (coordType == sphericalCoordSystem ? "Spherical" : "Unknown") ) )
               << " coordinate system." << std::endl;
    }

    // Get information about field variables

    // FIXME: Currently we assume that all fields live on the same
    // mesh. This may not be true.
    h5part_int64_t nFieldVars = H5BlockGetNumFields(file);
    if (nFieldVars < 0)
    {
        H5PartCloseFile(file);
        EXCEPTION2(NonCompliantFileException, "H5Part Constructor",
                   "Could not read number of field variables.");
    }
    else if (nFieldVars > 0)
    {
        for (h5part_int64_t idx=0; idx < nFieldVars; ++idx)
        {
            char varName[maxVarNameLen];
            h5part_int64_t gridRank;
            h5part_int64_t gridDims[3];
            h5part_int64_t fieldDims;
            h5part_int64_t type;

            if (H5BlockGetFieldInfo (file, idx, varName, maxVarNameLen,
                                     &gridRank, gridDims, &fieldDims, &type) !=
                H5PART_SUCCESS ) 
            {
                H5PartCloseFile(file);
                EXCEPTION2(NonCompliantFileException, "H5Part Constructor",
                           "Could not read field information.");
            }

            if (fieldDims == 1)
            {
                if (fieldScalarVarNameToTypeMap.find(varName) !=
                    fieldScalarVarNameToTypeMap.end())
                {
                    H5PartCloseFile(file);
                    EXCEPTION2(NonCompliantFileException, "H5Part Constructor",
                               "Field data contains two scalar variables with "
                               "the same name.");
                }
                else 
                {
                    fieldScalarVarNameToTypeMap[varName] = type;
                }
            }
            else
            {
                if (fieldVectorVarNameToTypeMap.find(varName) !=
                    fieldVectorVarNameToTypeMap.end())
                {
                    H5PartCloseFile(file);
                    EXCEPTION2(NonCompliantFileException, "H5Part Constructor",
                               "Field data contains two vector variables with "
                               "the same name.");
                }
                else 
                {
                    fieldVectorVarNameToTypeMap[varName] = type;
                    fieldVectorVarNameToFieldRankMap[varName] = fieldDims;
                }
            }
        }
    }

    // FIXME: This information should be read as attribute
    // from file.
    defaultIdVariableName = "id";
    
    h5part_int64_t numFileAttr = H5PartGetNumFileAttribs( file );
    
    for( h5part_int64_t i=0; i<numFileAttr; ++i )
    {
        char attrib_name[maxVarNameLen];
        h5part_int64_t len_of_attrib_name;
        h5part_int64_t attrib_type;
        h5part_int64_t attrib_nelem;

        H5PartGetFileAttribInfo ( file, i, attrib_name, maxVarNameLen,
                                  &attrib_type, &attrib_nelem );

        if( strncmp( attrib_name, "sortedKey", 9 ) == 0 )
        {
            // If the data was sorted use the sortedkey as the
            // default. Otherwise assume the id var name is "id".
            char sortedKey[maxVarNameLen];
            memset(sortedKey, 0, sizeof(char) * maxVarNameLen);
            
            if (H5PartReadFileAttrib(file, "sortedKey", sortedKey) == H5PART_SUCCESS)
            {
                if( std::string(sortedKey) != std::string("unsorted") )
                {
                    // FIXME: This information should be read as
                    // attribute from file.
                    defaultIdVariableName = sortedKey;
                }
            }
        }
    }

    idVariableName = defaultIdVariableName;
    
    // FIXME: Still need to check whether there are duplicate variable
    // names in field, scalar and vector variables and possibly modify
    // the name so that VisIt can distinguish between them.
    H5PartCloseFile(file);
    file = 0; // Mark file as closed
    
    visitTimer->StopTimer(t1, "H5PartFileFormat::avtH5PartFileFormat()");
}

// ****************************************************************************
//  Method: avtH5PartFileFormat destructor
//
//  Programmer: ghweber
//  Creation:   Tue eb 9 13:44:50 PST 2010
//
//  Modifications:
//    Gunther H. Weber, Wed Aug 17 16:27:07 PDT 2011
//    Close HDF5 file if this has not happened so far.
//
// ****************************************************************************

avtH5PartFileFormat::~avtH5PartFileFormat()
{
    FreeUpResources();
}

// ****************************************************************************
//  Method: avtEMSTDFileFormat::GetNTimesteps
//
//  Purpose:
//      Tells the rest of the code how many timesteps there are in this file.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
// Modifications:
//    Gunther H. Weber, Wed Aug  4 15:22:23 PDT 2010
//    Originally this function opened the file and used H5PartGetNumSteps()
//    to obtain the number of time steps. However, this call takes a very
//    long time and was slowing down VisIt considerably. To solve this
//    problem, we read the number of time steps when opening the file and
//    cache it as member variable.
//
// ****************************************************************************

int
avtH5PartFileFormat::GetNTimesteps(void)
{
    return numTimestepsInFile;;
}


// ****************************************************************************
//  Method: avtH5PartFileFormat::GetCycles
//
//  Purpose:
//      Returns the cycles
//
//  Arguments:
//      c          the cycles
//
//  Programmer: allen
//  Creation:   
//
// ****************************************************************************


void avtH5PartFileFormat::GetCycles(std::vector<int> &c)
{
    c = cycles;
}


// ****************************************************************************
//  Method: avtH5PartFileFormat::GetTimes
//
//  Purpose:
//      Returns the times
//
//  Arguments:
//      t          the times
//
//  Programmer: allen
//  Creation:   
//
// ****************************************************************************

void avtH5PartFileFormat::GetTimes(std::vector<double> &t)
{
    t = times;
}


// ****************************************************************************
//  Method: avtH5PartFileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it. This method is the mechanism for doing
//      that.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
// Modifications:
//    Gunther H. Weber, Wed Aug 17 13:14:27 PDT 2011
//    Honor useFastBitIndex option.
//
// ****************************************************************************

void
avtH5PartFileFormat::FreeUpResources(void)
{
    int t1 = visitTimer->StartTimer();

    if (file)
    {
        H5PartCloseFile(file);
        file = 0;
    }

    visitTimer->StopTimer(t1, "H5PartFileFormat::FreeUpResources()");
}


// ****************************************************************************
//  Method: avtH5PartFileFormat::DoubleToString
//
//  Purpose:
//      Create string representation of a double variable.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
// ****************************************************************************


std::string inline avtH5PartFileFormat::DoubleToString(double x)
{
    std::ostringstream o;
    if (!(o << setprecision(32) << x))
        EXCEPTION1(VisItException, "Error converting double to string.");
    return o.str();
}

// ****************************************************************************
//  Method: avtH5PartFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
// ****************************************************************************

void
avtH5PartFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md,
                                              int timeState)
{
    int t1 = visitTimer->StartTimer();

    if (particleVarNameToTypeMap.size())
    {
        // AddMeshToMetaData(md, meshname, mt, extents, nblocks, block_origin,
        //                   spatial_dimension, topological_dimension);
        std::string particleMeshName = "particles";
        AddMeshToMetaData(md, particleMeshName, AVT_POINT_MESH,
                          0 /* extents not known */,
                          1 /* single domain */,
                          0, particleNSpatialDims, 1);
        
        for (VarNameToInt64Map_t::const_iterator
               it = particleVarNameToTypeMap.begin();
             it != particleVarNameToTypeMap.end(); ++it)
        {
            // AddScalarVarToMetaData(md, varName, mesh_for_this_var, cent);
            AddScalarVarToMetaData(md, it->first, particleMeshName, AVT_NODECENT);
        }
    }

    // FIXME: Currently we assume that all fields live on the same
    // mesh. This may not be true.
    if (fieldScalarVarNameToTypeMap.size() ||
        fieldVectorVarNameToTypeMap.size())
    {
        std::string fieldMeshName = "fields";

        // AddMeshToMetaData(md, meshname, mt, extents, nblocks, block_origin,
        //                   spatial_dimension, topological_dimension);
        AddMeshToMetaData(md, fieldMeshName, AVT_RECTILINEAR_MESH,
                          0, 1, 0, 3, 3);
        // FIXME: Assume three dimensions for grid data
        // FIXME: Need to determine dimensions of mesh and create a
        // new mesh if dimensions for variables differ

        for (VarNameToInt64Map_t::const_iterator
               it = fieldScalarVarNameToTypeMap.begin();
             it != fieldScalarVarNameToTypeMap.end(); ++it)
        {
            // AddScalarVarToMetaData(md, varName, mesh_for_this_var, cent);
            AddScalarVarToMetaData(md, it->first, fieldMeshName, AVT_NODECENT);
        }

        for (VarNameToInt64Map_t::const_iterator
               it = fieldVectorVarNameToFieldRankMap.begin();
             it != fieldVectorVarNameToFieldRankMap.end(); ++it)
        {
            // AddVectorVarToMetaData(md, varName, mesh_for_this_var, cent,vector_dim);
            AddVectorVarToMetaData(md, it->first, fieldMeshName, AVT_NODECENT, it->second);
        }
    }

    if (enableDomainDecomposition)
        md->SetFormatCanDoDomainDecomposition(true);

    md->SetCyclesAreAccurate(true);
    md->SetCycles( cycles );

    if( times.size() )
    {
      md->SetTimesAreAccurate(true);
      md->SetTimes( times );
    }
    
    //
    // Here's the way to add expressions:
    //Expression momentum_expr;
    //momentum_expr.SetName("momentum");
    //momentum_expr.SetDefinition("{u, v}");
    //momentum_expr.SetType(Expression::VectorMeshVar);
    //md->AddExpression(&momentum_expr);
    //Expression KineticEnergy_expr;
    //KineticEnergy_expr.SetName("KineticEnergy");
    //KineticEnergy_expr.SetDefinition("0.5*(momentum*momentum)/(rho*rho)");
    //KineticEnergy_expr.SetType(Expression::ScalarMeshVar);
    //md->AddExpression(&KineticEnergy_expr);
    //
    visitTimer->StopTimer(t1, "H5PartFileFormat::PopulateDatabaseMetaData()");
}

// ****************************************************************************
//  Method: avtH5PartFileFormat::SelectParticlesToRead
//
//  Purpose:
//      Restrict the "view" in H5Part to those particles that we are
//      interested in: serial, no query -> all particles; serial, query ->
//      query result particles; parallel -> a subset of particles based on
//      query and domain decomposition settings.
//
//  Arguments:
//
//  Programmer: ghweber
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
//  Modifications:
//    Gunther H. Weber, Thu Apr  1 18:09:09 PDT 2010
//    Fixed computation of index range. Indices passed to H5PartSetView are
//    inclusive!
//
//    Gunther H. Weber, Wed Aug  4 19:18:23 PDT 2010
//    Move re-running of query to this method since it seems more logical
//    here.
//
//    Iuri Prilepov and Gunther H. Weber, Mon Aug 15 15:19:21 PDT 2011
//    Fix parallel decomposition of selections bug.
//
//    Eddie Rusu, Wed Oct 14 15:52:13 PDT 2020
//    Removed FastBit and FastQuery
//
// ****************************************************************************

void
avtH5PartFileFormat::SelectParticlesToRead( const char *varName )
{
    int t1 = visitTimer->StartTimer();

    // Reset the view so all particles are considered.
    H5PartSetView(file, -1, -1);

#ifdef PARALLEL
    h5part_int64_t nRanks = PAR_Size();
    h5part_int64_t rank   = PAR_Rank();
        
        // Read the whole file.
        h5part_int64_t nParticles = H5PartGetNumParticles(file);
        
        // Domain decomposition so read only the particles for this rank.
        if (enableDomainDecomposition)
        {
            // Determine what portion of the file will be read by each rank.
            h5part_int64_t idStart, idEnd, totalRead;
            GetDecomp( nParticles, idStart, idEnd, totalRead );

            // Restrict the view to the particles read on this rank.
            if( totalRead )
            {
                H5PartSetView(file, idStart, idEnd);
              
                debug5 << "SelectParticlesToRead(): Rank " << rank
                       << " reads " << totalRead << " of " << nParticles
                       << " particles from " << idStart << " to " << idEnd
                       << " in the file. " << std::endl;
            }
            // This rank does not have any particles to read.
            else
            {
                h5part_int64_t *indices = new h5part_int64_t[0];
                H5PartSetViewIndices(file, indices, 0);
                delete[] indices;
                
                debug5 << "SelectParticlesToRead(): Rank " << rank
                       << " no particles to read." << std::endl;
            }            
        }
        // Read the entire data set on rank zero
        else
        {
            debug5 << "SelectParticlesToRead(): Rank " << rank
                   << " domain decomposition disabled." << std::endl;
            
            if( PAR_Rank() )
            {
                h5part_int64_t *indices = new h5part_int64_t[0];
                H5PartSetViewIndices(file, indices, 0);
                delete[] indices;

                debug5 << "SelectParticlesToRead(): Rank " << rank
                       << " no particles to read." << std::endl;
            }
            else
            {
                H5PartSetView(file, -1, -1);
          
                debug5 << "SelectParticlesToRead(): Rank " << rank
                       << " All " << nParticles << " particles "
                       << " in the file read." << std::endl;
            }
        }

#else
        // Read the whole file.
        H5PartSetView(file, -1, -1);

        h5part_int64_t nParticles = H5PartGetNumParticles(file);

        debug5 << "SelectParticlesToRead(): Serial "
               << " reading all particles (" << nParticles
               << ") in file (no query)." << std::endl;
#endif

    visitTimer->StopTimer(t1, "H5PartFileFormat::SelectParticlesToRead()");
}


// ****************************************************************************
//  Method: avtH5PartFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      timestate   The index of the timestate.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
//  Modifications:
//
// ****************************************************************************

vtkDataSet *
avtH5PartFileFormat::GetMesh(int timestate, const char *meshname)
{
    if (strcmp(meshname, "particles") == 0)
    {
        return GetParticleMesh(timestate); // There is only one particle mesh
    }
    else
    {
        return GetFieldMesh(timestate, meshname); // In the future there may be
        // more than one mesh for field data
    }
}

// ****************************************************************************
//  Method: avtH5PartFileFormat::GetParticleMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      timestate   The index of the timestate.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
//  Modifications:
//    Brad Whitlock, Mon Dec 12 16:06:13 PST 2011
//    Construct avtOriginalCellNumbers here if we're doing a data selection.
//
// ****************************************************************************

vtkDataSet *
avtH5PartFileFormat::GetParticleMesh(int timestate)
{
    int t1 = visitTimer->StartTimer();

    // Switch to appropriate time step (opens file if necessary)
    ActivateTimestep(timestate);

    idVariableName = defaultIdVariableName;

    VarNameToInt64Map_t::const_iterator it =
      particleVarNameToTypeMap.find(idVariableName);
    
    if (it == particleVarNameToTypeMap.end())
      EXCEPTION1(InvalidVariableException, idVariableName);
    
    // Select particles to actually read
    SelectParticlesToRead( idVariableName.c_str() );
    
    // Read data
    h5part_int64_t nPoints = H5PartGetNumParticles(file);

    // If we do domain decomposition or particle selection, this will only
    // return the number of particles in the current view, i.e., the number of
    // particles we actually read
    debug5 << "GetParticleMesh(): Reading " << nPoints
           << " particles." << std::endl;

    if (nPoints ==  0)
    {
        visitTimer->StopTimer(t1, "H5PartFileFormat::GetParticleMesh()");
        return 0;
    }

    const char *coordNames[3][3] = {
        { "x", "y", "z" }, {"r", "phi", "z" } , { "r", "phi", "theta" }
    };

    h5part_int64_t coordValType = particleVarNameToTypeMap[coordNames[coordType][0]]; 
    if (coordValType !=
        (h5part_int64_t) particleVarNameToTypeMap[coordNames[coordType][1]] ||
        ((particleNSpatialDims > 2) &&
         coordValType != (h5part_int64_t) particleVarNameToTypeMap[coordNames[coordType][2]])) /// TODO: check parenthsis placement logic 
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetParticleMesh",
                   "Coordinate data sets do not have the same value type "
                   "(double, float, ...).");
    }

    void* coordsRead[3] = { 0, 0, 0 };

    for (int i=0; i < particleNSpatialDims; ++i)
    {
        h5part_int64_t status = H5PART_SUCCESS; 

        if (coordValType == H5PART_INT64)
        {
            coordsRead[i] = malloc(sizeof(h5part_int64_t) * nPoints);
            status = H5PartReadDataInt64(file, coordNames[coordType][i],
                    static_cast<h5part_int64_t*>(coordsRead[i]));
        }
        else if (coordValType == H5PART_INT32)
        {
            coordsRead[i] = malloc(sizeof(h5part_int32_t) * nPoints);
            status = H5PartReadDataInt32(file, coordNames[coordType][i],
                    static_cast<h5part_int32_t*>(coordsRead[i]));
        }
        else if (coordValType == H5PART_FLOAT64)
        {
            coordsRead[i] = malloc(sizeof(h5part_float64_t) * nPoints);
            status = H5PartReadDataFloat64(file, coordNames[coordType][i],
                    static_cast<h5part_float64_t*>(coordsRead[i]));

        }
        else if (coordValType == H5PART_FLOAT32)
        {
            coordsRead[i] = malloc(sizeof(h5part_float32_t) * nPoints);
            status = H5PartReadDataFloat32(file, coordNames[coordType][i],
                    static_cast<h5part_float32_t*>(coordsRead[i]));
        }
        else
        {
            EXCEPTION2(NonCompliantFileException, "H5Part GetParticleMesh",
                       "Unsupported value type for coordinates. "
                       "(Supported tyes are FLOAT64, FLOAT32 and INT64.)");
        }
        // FIXME: Possibly handle other data types for coordinates.

        if (status != H5PART_SUCCESS)
        {
            EXCEPTION2(NonCompliantFileException, "H5Part GetParticleMesh",
                       "Could not read coordinates.");
        }
    }

    // Construct VTK data set
    vtkUnstructuredGrid *dataset = vtkUnstructuredGrid::New();
    vtkPoints *vtkpoints = vtkPoints::New();

    // FIXME: Once VisIt supports dobule data, we need to ensure that
    // native precsission is preserved
    vtkpoints->SetNumberOfPoints((vtkIdType) nPoints);

    float *pts = (float *) vtkpoints->GetVoidPointer(0);

    if (coordValType == H5PART_INT64)
    {
        for (h5part_int64_t ptNo = 0; ptNo < nPoints; ++ptNo)
        {
            pts[3*ptNo] =
                float(static_cast<h5part_int64_t*>(coordsRead[0])[ptNo]);
            pts[3*ptNo+1] =
                float(static_cast<h5part_int64_t*>(coordsRead[1])[ptNo]);
            if (particleNSpatialDims > 2)
                pts[3*ptNo+2] =
                    float(static_cast<h5part_int64_t*>(coordsRead[2])[ptNo]);
            else
                pts[3*ptNo+2] = 0.0;
        }
    }
    else if (coordValType == H5PART_INT32)
    {
        for (h5part_int32_t ptNo = 0; ptNo < nPoints; ++ptNo)
        {
            pts[3*ptNo] =
                float(static_cast<h5part_int32_t*>(coordsRead[0])[ptNo]);
            pts[3*ptNo+1] =
                float(static_cast<h5part_int32_t*>(coordsRead[1])[ptNo]);
            if (particleNSpatialDims > 2)
                pts[3*ptNo+2] =
                    float(static_cast<h5part_int32_t*>(coordsRead[2])[ptNo]);
            else
                pts[3*ptNo+2] = 0.0;
        }
    }
    else if (coordValType == H5PART_FLOAT64)
    {
        for (h5part_int64_t ptNo = 0; ptNo < nPoints; ++ptNo)
        {
            pts[3*ptNo] =
                float(static_cast<h5part_float64_t*>(coordsRead[0])[ptNo]);
            pts[3*ptNo+1] =
                float(static_cast<h5part_float64_t*>(coordsRead[1])[ptNo]);
            if (particleNSpatialDims > 2)
                pts[3*ptNo+2] =
                    float(static_cast<h5part_float64_t*>(coordsRead[2])[ptNo]);
            else
                pts[3*ptNo+2] = 0.0;
        }
    }
    else if (coordValType == H5PART_FLOAT32)
    {
        for (h5part_int64_t ptNo = 0; ptNo < nPoints; ptNo++)
        {
            pts[3*ptNo] =
                float(static_cast<h5part_float32_t*>(coordsRead[0])[ptNo]);
            pts[3*ptNo+1] =
                float(static_cast<h5part_float32_t*>(coordsRead[1])[ptNo]);
            if (particleNSpatialDims > 2)
                pts[3*ptNo+2] =
                    float(static_cast<h5part_float32_t*>(coordsRead[2])[ptNo]);
            else
                pts[3*ptNo+2] = 0.0;
        }
    }
    else
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetParticleMesh",
                   "Unsupported value type for coordinates. "
                   "Supported types are FLOAT64. FLOAT32 and INT64.");
    }
    // FIXME: Possibly handle other data types for the coordinates

    for (int i=0; i <particleNSpatialDims; ++i)
        free(coordsRead[i]);

    dataset->Allocate(nPoints);

    for (h5part_int64_t i=0; i<nPoints; ++i)
    {
        vtkIdType onevertex = (vtkIdType) i;
        dataset->InsertNextCell(VTK_VERTEX, 1, &onevertex);
    }

    dataset->SetPoints(vtkpoints);
    vtkpoints->Delete();

    //
    // If a data selection was processed successfully then dataset
    // with different connectivity from the original full dataset will
    // be produced.  This means that if we're requesting original cell
    // numbers above up in the database then the values calculated
    // there would be wrong. We need to calculate them here and send
    // them along so named selections will work properly.
    //
    debug5 << "avtH5PartFileFormat::GetParticleMesh(): "
           << " Creating avtOriginalCellNumbers early" << std::endl;

    vtkUnsignedIntArray *origZones = vtkUnsignedIntArray::New();
    origZones->SetName("avtOriginalCellNumbers");
    origZones->SetNumberOfComponents(2);
    origZones->SetNumberOfTuples(nPoints);
    unsigned int *iptr = (unsigned int *) origZones->GetVoidPointer(0);

    if (it->second == H5PART_FLOAT64)
    {
        h5part_float64_t *idList = new h5part_float64_t[nPoints];;
        if (H5PartReadDataFloat64(file, idVariableName.c_str(), idList) != H5PART_SUCCESS)
        {
          delete[] idList;
          EXCEPTION1(InvalidVariableException, idVariableName);
        }
        else
        {
          for(size_t i = 0; i < nPoints; ++i)
          {
            // domain 0. This format is single domain so it's probably okay.
            *iptr++ = 0;
            *iptr++ = (unsigned int) idList[i];
        
            debug5 << "Domain " << 0 << " cell " << idList[i]
                   << std::endl;
          }
    
          delete[] idList;
        }
    }
    else if (it->second == H5PART_FLOAT32)
    {
        h5part_float32_t *idList = new h5part_float32_t[nPoints];
        if (H5PartReadDataFloat32(file, idVariableName.c_str(), idList) != H5PART_SUCCESS)
        {
          delete[] idList;
          EXCEPTION1(InvalidVariableException, idVariableName);
        }
        else
        {
          for(size_t i = 0; i < nPoints; ++i)
          {
            // domain 0. This format is single domain so it's probably okay.
            *iptr++ = 0;
            *iptr++ = (unsigned int) idList[i];
        
            debug5 << "Domain " << 0 << " cell " << idList[i]
                   << std::endl;
          }

          delete[] idList;
        }
    }
    else if (it->second == H5PART_INT64)
    {
      h5part_int64_t *idList = new h5part_int64_t[nPoints];
      if (H5PartReadDataInt64(file, idVariableName.c_str(), idList) != H5PART_SUCCESS)
      {
          delete[] idList;
          EXCEPTION1(InvalidVariableException, idVariableName);
      }
      else
      {
          for(size_t i = 0; i < nPoints; ++i)
          {
            // domain 0. This format is single domain so it's probably okay.
            *iptr++ = 0;
            *iptr++ = (unsigned int) idList[i];
        
            debug5 << "Domain " << 0 << " cell " << idList[i]
                   << std::endl;
          }

          delete[] idList;
      }
    }
    else if (it->second == H5PART_INT32)
    {
      h5part_int32_t *idList = new h5part_int32_t[nPoints];
      if (H5PartReadDataInt32(file, idVariableName.c_str(), idList) != H5PART_SUCCESS)
      {
          delete[] idList;
          EXCEPTION1(InvalidVariableException, idVariableName);
      }
      else
      {
          for(size_t i = 0; i < nPoints; ++i)
          {
            // domain 0. This format is single domain so it's probably okay.
            *iptr++ = 0;
            *iptr++ = (unsigned int) idList[i];
        
            debug5 << "Domain " << 0 << " cell " << idList[i]
                   << std::endl;
          }

          delete[] idList;
      }
    }
    else
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetVar",
                   std::string("Scalar variable ") + idVariableName + 
                   std::string(" uses an unsupported value type. Supported "
                               "types are FLOAT64, FLOAT32 and INT64."));
    }
    // FIXME: Possibly handle other data types for the id variable

    dataset->GetCellData()->AddArray(origZones);
    origZones->Delete();
    
    visitTimer->StopTimer(t1, "H5PartFileFormat::GetParticleMesh()");

    // Reset view
    H5PartSetView(file, -1, -1);

    return dataset;
}

// ****************************************************************************
//  Method: avtH5PartFileFormat::GetFieldMesh
//
//  Purpose:
//      Gets the mesh associated with fields of this file
//
//  Arguments:
//      timestate   The index of the timestate.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      domain      The index of the domain.  If there are NDomains, this
//                  value is guaranteed to be between 0 and NDomains-1,
//                  regardless of block origin.
//
//  Programmer: kurts
//  Creation:   Tue Aug 28 17:35:50 PDT 2007
//
//  Modifications:
//    Gunther H. Weber, Fri Apr 17 13:55:07 PDT 2009
//    Read block origin from file.
//
//    Gunther H. Weber, Fri Jul 17 11:11:36 PDT 2009
//    Clean up. Read block origin via dedicated API call.
//
//    Gunther H. Weber, Wed Feb 10 11:47:32 PST 2010
//    Conversion from MD plugin to SD plugin supporting domain decomposition
//    and clean-up.
//
//    Gunther H. Weber, Thu Apr  1 18:14:00 PDT 2010
//    Adapted to GetSubBlock paramter change.
//
// ****************************************************************************

vtkDataSet *
avtH5PartFileFormat::GetFieldMesh(int timestate, const char *meshname)
{
    int t1 = visitTimer->StartTimer();

    // Switch to appropriate time step (opens file if necessary)
    ActivateTimestep(timestate);

    // FIXME: assume all fields have the same dimensions
    // get information from first field
    char fieldName[maxVarNameLen];
    h5part_int64_t gridRank;
    h5part_int64_t gridDims[3];
    h5part_int64_t fieldDims;
    h5part_int64_t type;
    if (H5BlockGetFieldInfo (file, 0, fieldName, maxVarNameLen,
                &gridRank, gridDims, &fieldDims, &type) != H5PART_SUCCESS)
        EXCEPTION2(NonCompliantFileException, "H5Part GetFieldMesh",
                   "Could not read field information.");

    // Determine which data portion we need to construct a mesh for
    h5part_int64_t subBlockDims[6];
    GetSubBlock(gridDims, subBlockDims);

    double xOrigin = 0;
    double yOrigin = 0;
    double zOrigin = 0;

    if (H5Block3dGetFieldOrigin(file, fieldName,
                &xOrigin, &yOrigin, &zOrigin) != H5PART_SUCCESS)
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetFieldMesh",
                   "Could not read field origin.");
    }

    debug5 << "xOrigin: " << xOrigin << "\n";
    debug5 << "yOrigin: " << yOrigin << "\n";
    debug5 << "zOrigin: " << zOrigin << "\n";

    double xSpacing = 0;
    double ySpacing = 0;
    double zSpacing = 0;

    if (H5Block3dGetFieldSpacing (file, fieldName,
                &xSpacing, &ySpacing, &zSpacing) != H5PART_SUCCESS)
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetFieldMesh", 
                   "Could not read field spacing.");
    }

    debug5 << "xSpacing: " << xSpacing << "\n";
    debug5 << "ySpacing: " << ySpacing << "\n";
    debug5 << "zSpacing: " << zSpacing << "\n";

    vtkFloatArray *coords[3] = {0, 0, 0};
    // set x coordinates
    coords[0] = vtkFloatArray::New();
    coords[0]->SetNumberOfTuples(subBlockDims[1] - subBlockDims[0] + 1);
    float *xarray = (float *) coords[0]->GetVoidPointer(0);
    for (int i=subBlockDims[0]; i <= subBlockDims[1]; i++)
    {
        xarray[i-subBlockDims[0]] = xOrigin + i * xSpacing;
    }

    // set y coordinates
    coords[1] = vtkFloatArray::New();
    coords[1]->SetNumberOfTuples(subBlockDims[3] - subBlockDims[2] + 1);
    float *yarray = (float *) coords[1]->GetVoidPointer(0);
    for (int i=subBlockDims[2]; i <= subBlockDims[3]; i++)
    {
        yarray[i-subBlockDims[2]] = yOrigin + i * ySpacing;
    }

    // set z coordinates
    coords[2] = vtkFloatArray::New();
    coords[2]->SetNumberOfTuples(subBlockDims[5] - subBlockDims[4] + 1);
    float *zarray = (float *) coords[2]->GetVoidPointer(0);
    for (int i=subBlockDims[4]; i <= subBlockDims[5]; i++)
    {
        zarray[i-subBlockDims[4]] = zOrigin + i * zSpacing;
    }

    // create vtkRectilinearGrid objects + set dims and coords
    vtkRectilinearGrid *rgrid = vtkRectilinearGrid::New();
    rgrid->SetDimensions(
            coords[0]->GetNumberOfTuples(),
            coords[1]->GetNumberOfTuples(),
            coords[2]->GetNumberOfTuples());

    rgrid->SetXCoordinates(coords[0]);
    coords[0]->Delete();
    rgrid->SetYCoordinates(coords[1]);
    coords[1]->Delete();
    rgrid->SetZCoordinates(coords[2]);
    coords[2]->Delete();

    visitTimer->StopTimer(t1, "H5PartFileFormat::GetFieldMesh()");

    return rgrid;
}

// ****************************************************************************
//  Method: avtH5PartFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      varName    The name of the variable requested.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
// ****************************************************************************

vtkDataArray *
avtH5PartFileFormat::GetVar(int timestate, const char *varName)
{
    int t1 = visitTimer->StartTimer();

    VarNameToInt64Map_t::const_iterator it =
      particleVarNameToTypeMap.find(varName);

    if (it == particleVarNameToTypeMap.end())
    {
        visitTimer->StopTimer(t1, "H5PartFileFormat::GetVar() [invoking GetFieldVar]");
        return GetFieldVar(timestate, varName);
    }

    // Switch to appropriate time step (opens file if necessary)
    ActivateTimestep(timestate);

    // Select particles to actually read
    SelectParticlesToRead( varName );

    // Read data
    h5part_int64_t nPoints = H5PartGetNumParticles(file);

    debug5 << "GetVar() reading " << nPoints << " particles." << std::endl;

    if (nPoints == 0)
    {
        visitTimer->StopTimer(t1, "H5PartFileFormat::GetVar()");
        return 0;
    }

    vtkDataArray *scalars = 0;
    if (it->second == H5PART_FLOAT64)
    {
        scalars = vtkDoubleArray::New();
        scalars->SetNumberOfTuples(nPoints);
        if (H5PartReadDataFloat64(file, varName,
                    static_cast<h5part_float64_t*>(scalars->GetVoidPointer(0)))
                != H5PART_SUCCESS)
        {
            EXCEPTION1(InvalidVariableException, varName);
        }
    }
    else if (it->second == H5PART_FLOAT32)
    {
        scalars = vtkFloatArray::New();
        scalars->SetNumberOfTuples(nPoints);
        if (H5PartReadDataFloat32(file, varName,
                    static_cast<h5part_float32_t*>(scalars->GetVoidPointer(0)))
                != H5PART_SUCCESS)
        {
            EXCEPTION1(InvalidVariableException, varName);
        }
    }
    else if (it->second == H5PART_INT64)
    {
        // Unlike FLOAT32 and FLOAT64 which should map to float and
        // double we are not sure what int type is used to achieve 64 bits.
        // To be safe, read into separate array and copy.
        h5part_int64_t *idvar = new h5part_int64_t[nPoints];
        if (H5PartReadDataInt64(file, varName, idvar) != H5PART_SUCCESS)
            EXCEPTION1(InvalidVariableException, varName);

        scalars = vtkFloatArray::New(); // FIXME: Can an integer array be used?
        scalars->SetNumberOfTuples(nPoints);
        float *ptr = static_cast<float*>(scalars->GetVoidPointer(0));
        for (h5part_int64_t i=0; i<nPoints; ++i)
        {
            ptr[i] = float(idvar[i]);
        }
        delete[] idvar;
    }
    else if (it->second == H5PART_INT32)
    {
        // Unlike FLOAT32 and FLOAT64 which should map to float and
        // double we are not sure what int type is used to achieve 32 bits.
        // To be safe, read into separate array and copy.
        h5part_int32_t *idvar = new h5part_int32_t[nPoints];
        if (H5PartReadDataInt32(file, varName, idvar) != H5PART_SUCCESS)
            EXCEPTION1(InvalidVariableException, varName);

        scalars = vtkFloatArray::New(); // FIXME: Can an integer array be used?
        scalars->SetNumberOfTuples(nPoints);
        float *ptr = static_cast<float*>(scalars->GetVoidPointer(0));
        for (h5part_int32_t i=0; i<nPoints; ++i)
        {
            ptr[i] = float(idvar[i]);
        }
        delete[] idvar;
    }
    else
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetVar",
                   std::string("Scalar variable ") + varName + 
                   std::string(" uses an unsupported value type. Supported "
                               "types are FLOAT64, FLOAT32 and INT64."));
    }
    // FIXME: Possibly handle other data types for scalar variables

    // Reset view
    H5PartSetView(file, -1, -1);

    visitTimer->StopTimer(t1, "H5PartFileFormat::GetVar()");

    return scalars;
}


// ****************************************************************************
//  Method: avtH5PartFileFormat::GetFieldVar
//
//  Purpose:
//      Get a scalar variable associated with the field mesh of this file
//
//  Arguments:
//      timestate   The index of the timestate.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      varName    The name of the variable requested.
//
//  Programmer: kurts
//  Creation:   Tue Aug 28 17:35:50 PDT 2007
//
//  Modifications:
//    Gunther H. Weber, Fri Jul 17 11:13:39 PDT 2009
//    Clean-up.
//
//    Gunther H. Weber, Wed Feb 10 11:47:32 PST 2010
//    Conversion from MD plugin to SD plugin supporting domain decomposition.
//
//    Gunther H. Weber, Thu Apr  1 18:14:00 PDT 2010
//    Adapted to GetSubBlock paramter change.
//
// ****************************************************************************

vtkDataArray *
avtH5PartFileFormat::GetFieldVar(int timestate, const char* varName)
{
    int t1 = visitTimer->StartTimer();

    // Check if we know about variable and get type
    VarNameToInt64Map_t::const_iterator it =
      fieldScalarVarNameToTypeMap.find(varName);

    if (it == fieldScalarVarNameToTypeMap.end())
        EXCEPTION1(InvalidVariableException, varName);

    // Activate correct time step (opens file if necessary)
    ActivateTimestep(timestate);

    char fieldName[maxVarNameLen];
    h5part_int64_t gridRank;
    h5part_int64_t gridDims[3];
    h5part_int64_t fieldDims;
    h5part_int64_t type;

    // FIXME: This assumes that all blocks have the same mesh
    if (H5BlockGetFieldInfo (file, 0, fieldName, maxVarNameLen,
                &gridRank, gridDims, &fieldDims, &type) != H5PART_SUCCESS)
        EXCEPTION1(InvalidVariableException,
                   "Could not read field information.");

    h5part_int64_t subBlockDims[6];
    GetSubBlock(gridDims, subBlockDims);

    // set field layout
    if(H5BlockDefine3DFieldLayout(file, subBlockDims[0], subBlockDims[1], 
                subBlockDims[2], subBlockDims[3], 
                subBlockDims[4], subBlockDims[5]) != H5PART_SUCCESS)
        EXCEPTION1(VisItException, "Could not set field layout.");

    // read data
    int nValues = (subBlockDims[1] - subBlockDims[0] + 1) *
        (subBlockDims[3] - subBlockDims[2] + 1) *
        (subBlockDims[5] - subBlockDims[4] + 1);

    vtkDataArray *scalars;
    if (it->second == H5PART_FLOAT64)
    {
        scalars = vtkDoubleArray::New();
        scalars->SetNumberOfTuples(nValues);
        if (H5Block3dReadScalarFieldFloat64(file, varName,
                    static_cast<h5part_float64_t*>(scalars->GetVoidPointer(0)))
                != H5PART_SUCCESS)
        {
            EXCEPTION1(InvalidVariableException, varName);
        }
    }
    else if (it->second == H5PART_FLOAT32)
    {
        scalars = vtkFloatArray::New();
        scalars->SetNumberOfTuples(nValues);
        if (H5Block3dReadScalarFieldFloat32(file, varName,
                    static_cast<h5part_float32_t*>(scalars->GetVoidPointer(0)))
                != H5PART_SUCCESS)
        {
            EXCEPTION1(InvalidVariableException, varName);
        }
    }
    else if (it->second == H5PART_INT64)
    {
        // Unlike FLOAT32 and FLOAT64 which should map to float and
        // double we are not sure what int type is used to achieve 64 bits.
        // To be safe, read into separate array and copy.
        h5part_int64_t *idvar = new h5part_int64_t[nValues];
        if (H5Block3dReadScalarFieldInt64(file, varName, idvar) != H5PART_SUCCESS)
            EXCEPTION1(InvalidVariableException, varName);

        scalars = vtkFloatArray::New(); // FIXME: Can an integer array be used?
        scalars->SetNumberOfTuples(nValues);
        float *ptr = static_cast<float*>(scalars->GetVoidPointer(0));
        for (h5part_int64_t i=0; i<nValues; ++i)
        {
            ptr[i] = float(idvar[i]);
        }
        delete[] idvar;
    }
    else if (it->second == H5PART_INT32)
    {
        // Unlike FLOAT32 and FLOAT64 which should map to float and
        // double we are not sure what int type is used to achieve 32 bits.
        // To be safe, read into separate array and copy.
        h5part_int32_t *idvar = new h5part_int32_t[nValues];
        if (H5Block3dReadScalarFieldInt32(file, varName, idvar) != H5PART_SUCCESS)
            EXCEPTION1(InvalidVariableException, varName);

        scalars = vtkFloatArray::New(); // FIXME: Can an integer array be used?
        scalars->SetNumberOfTuples(nValues);
        float *ptr = static_cast<float*>(scalars->GetVoidPointer(0));
        for (h5part_int32_t i=0; i<nValues; ++i)
        {
            ptr[i] = float(idvar[i]);
        }
        delete[] idvar;
    }
    else
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetVar",
                   std::string("Field variable ") + varName + 
                   std::string(" uses an unsupported value type. Supported "
                               "types are FLOAT64, FLOAT32 and INT64."));
    }
    // FIXME: Possibly handle other data types for field variables

    visitTimer->StopTimer(t1, "H5PartFileFormat::GetFieldVar()");

    return scalars;
}

// ****************************************************************************
//  Method: avtH5PartFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      varName    The name of the variable requested.
//
//  Programmer: ghweber -- generated by xml2avt
//  Creation:   Tue Feb 9 13:44:50 PST 2010
//
//  Modifications:
//    Kurt Stockinger, Tue Aug 28 17:35:50 PDT 2007
//    Added support for field data
//
//    Gunther H. Weber, Fri Jul 17 11:13:39 PDT 2009
//    Clean-up.
//
//    Gunther H. Weber, Thu Apr  1 18:14:00 PDT 2010
//    Adapted to GetSubBlock paramter change.
//
// ****************************************************************************
vtkDataArray *
avtH5PartFileFormat::GetVectorVar(int timestate, const char *varName)
{
    int t1 = visitTimer->StartTimer();

    // Check if we know about variable and get type
    VarNameToInt64Map_t::const_iterator it =
      fieldVectorVarNameToTypeMap.find(varName);
    
    if (it == fieldVectorVarNameToTypeMap.end())
        EXCEPTION1(InvalidVariableException, varName);

    // Activate correct time step (opens file if necessary)
    ActivateTimestep(timestate);

    char fieldName[maxVarNameLen];
    h5part_int64_t gridRank;
    h5part_int64_t gridDims[3];
    h5part_int64_t fieldDims;
    h5part_int64_t type;

    // FIXME: This assumes that all blocks have the same size
    if (H5BlockGetFieldInfo (file, 0, fieldName, maxVarNameLen,
                &gridRank, gridDims, &fieldDims, &type) != H5PART_SUCCESS)
        EXCEPTION1(InvalidVariableException, "Could not read field information.");

    h5part_int64_t subBlockDims[6];
    GetSubBlock(gridDims, subBlockDims);

    // set field layout
    if (H5BlockDefine3DFieldLayout (file, 
                subBlockDims[0], subBlockDims[1], 
                subBlockDims[2], subBlockDims[3], 
                subBlockDims[4], subBlockDims[5]) != H5PART_SUCCESS)
        EXCEPTION1(InvalidVariableException, "Could not set field layout");

    // Read data
    int nValues = (subBlockDims[1] - subBlockDims[0] + 1) *
        (subBlockDims[3] - subBlockDims[2] + 1) *
        (subBlockDims[5] - subBlockDims[4] + 1);

    void* vecsRead[3] = { 0, 0, 0 };
    h5part_int64_t status = H5PART_SUCCESS; (void) status;
    if (it->second == H5PART_INT64)
    {
        for (int i=0; i<3; ++i)
            vecsRead[i] = malloc(sizeof(h5part_int64_t) * nValues);

        status = H5Block3dRead3dVectorFieldInt64(file, varName,  
                static_cast<h5part_int64_t*>(vecsRead[0]),
                static_cast<h5part_int64_t*>(vecsRead[1]),
                static_cast<h5part_int64_t*>(vecsRead[2]));
    }
    else if (it->second == H5PART_INT32)
    {
        for (int i=0; i<3; ++i)
            vecsRead[i] = malloc(sizeof(h5part_int32_t) * nValues);

        status = H5Block3dRead3dVectorFieldInt32(file, varName,  
                static_cast<h5part_int32_t*>(vecsRead[0]),
                static_cast<h5part_int32_t*>(vecsRead[1]),
                static_cast<h5part_int32_t*>(vecsRead[2]));
    }
    else if (it->second == H5PART_FLOAT64)
    {
        for (int i=0; i<3; ++i)
            vecsRead[i] = malloc(sizeof(h5part_float64_t) * nValues);

        status = H5Block3dRead3dVectorFieldFloat64(file, varName,  
                static_cast<h5part_float64_t*>(vecsRead[0]),
                static_cast<h5part_float64_t*>(vecsRead[1]),
                static_cast<h5part_float64_t*>(vecsRead[2]));


    }
    else if (it->second == H5PART_FLOAT32)
    {
        for (int i=0; i<3; ++i)
            vecsRead[i] = malloc(sizeof(h5part_float32_t) * nValues);

        status = H5Block3dRead3dVectorFieldFloat32(file, varName,  
                static_cast<h5part_float32_t*>(vecsRead[0]),
                static_cast<h5part_float32_t*>(vecsRead[1]),
                static_cast<h5part_float32_t*>(vecsRead[2]));
    }
    else
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetVectorVar",
                   std::string("Vector variable ") + varName +
                   std::string(" uses an unsupported value type. Supported"
                               "types are FLOAT64, FLOAT32 and INT64."));
    }
    // FIXME: Possibly handle other data types for vector variables

    // Store data in VTK object
    vtkDataArray *array;
    if (it->second == H5PART_FLOAT64)
    {
        array = vtkDoubleArray::New();
    }
    else
    {
        array = vtkFloatArray::New(); //TODO: check fix
        // FIXME: Possibly check for integer and create an integer array
    }

    // Set array size
    array->SetNumberOfComponents(3);
    array->SetNumberOfTuples(nValues);

    // Copy data
    if (it->second == H5PART_INT64)
    {
        // FIXME: Possibly check for integer and create an integer array
        float *vecs = static_cast<float*>(array->GetVoidPointer(0));
        for (h5part_int64_t valNo = 0; valNo < nValues; valNo++)
        {
            vecs[3*valNo] =
                float(static_cast<h5part_int64_t*>(vecsRead[0])[valNo]);
            vecs[3*valNo+1] =
                float(static_cast<h5part_int64_t*>(vecsRead[1])[valNo]);
            vecs[3*valNo+2] =
                float(static_cast<h5part_int64_t*>(vecsRead[2])[valNo]);
        }
    }
    else if (it->second == H5PART_INT32)
    {
        // FIXME: Possibly check for integer and create an integer array
        float *vecs = static_cast<float*>(array->GetVoidPointer(0));
        for (h5part_int32_t valNo = 0; valNo < nValues; valNo++)
        {
            vecs[3*valNo] =
                float(static_cast<h5part_int32_t*>(vecsRead[0])[valNo]);
            vecs[3*valNo+1] =
                float(static_cast<h5part_int32_t*>(vecsRead[1])[valNo]);
            vecs[3*valNo+2] =
                float(static_cast<h5part_int32_t*>(vecsRead[2])[valNo]);
        }
    }
    else if (it->second == H5PART_FLOAT64)
    {
        double *vecs = static_cast<double*>(array->GetVoidPointer(0));
        for (h5part_int64_t valNo = 0; valNo < nValues; valNo++)
        {
            vecs[3*valNo] =
                double(static_cast<h5part_float64_t*>(vecsRead[0])[valNo]);
            vecs[3*valNo+1] =
                double(static_cast<h5part_float64_t*>(vecsRead[1])[valNo]);
            vecs[3*valNo+2] =
                double(static_cast<h5part_float64_t*>(vecsRead[2])[valNo]);
        }
    }
    else if (it->second == H5PART_FLOAT32)
    {
        float *vecs = static_cast<float*>(array->GetVoidPointer(0));
        for (h5part_int64_t valNo = 0; valNo < nValues; valNo++)
        {
            vecs[3*valNo] =
                float(static_cast<h5part_float32_t*>(vecsRead[0])[valNo]);
            vecs[3*valNo+1] =
                float(static_cast<h5part_float32_t*>(vecsRead[1])[valNo]);
            vecs[3*valNo+2] =
                float(static_cast<h5part_float32_t*>(vecsRead[2])[valNo]);
        }
    }
    else
    {
        EXCEPTION2(NonCompliantFileException, "H5Part GetVectorVar",
                   std::string("Vector variable ") + varName +
                   std::string(" uses an unsupported value type. Supported "
                               "types are FLOAT64, FLOAT32 and INT64.)"));
    }
    // FIXME: Possibly handle other data types for vector variables
 
    for (int i=0; i<3; ++i)
      free(vecsRead[i]);

    visitTimer->StopTimer(t1, "H5PartFileFormat::GetVectorVar()");
    return array;
}

// ****************************************************************************
//  Method: avtH5PartFileFormat::ActivateTimestep
//
//  Purpose:
//      Select the current time step
//
//  Arguments:
//      ts      The new current time step
//
//  Output:
//
//  Programmer: ghweber
//  Creation:   Thu Feb 11 15:38:49 PST 2010
//
//  Modifications:
//    Gunther H. Weber, Mon Aug  9 14:58:58 PDT 2010
//    Fixed bug identified by regression tests that sometimes (it seems
//    mainly for field data) H5PartSetStep is not called when a new file
//    is opened.
//
//    Gunther H. Weber, Wed Aug 17 13:44:13 PDT 2011
//    Do not open file with HDF5_FQ if useFastBitIndex is false.
//
//    Eddie Rusu, Wed Oct 14 15:52:13 PDT 2020
//    Removed FastBit and FastQuery
//
// ****************************************************************************

void
avtH5PartFileFormat::ActivateTimestep(int ts)
{
    int t1 = visitTimer->StartTimer();

    // Open file if necessary
    if (!file)
    {
        debug5 << "avtH5PartFileFormat::ActivateTimestep(): Opening file and ";
        debug5 << "activating time step " << ts << std::endl;

        file = H5PartOpenFile(filenames[0], H5PART_READ);

        if (!file)
          EXCEPTION1(InvalidFilesException, "Cannot open file.");
        
        if (H5PartSetStep(file, ts) != H5PART_SUCCESS)
        {
          std::ostringstream msg;
          msg << "Cannot activate time step " << ts << ".";

          debug1 << "avtH5PartFileFormat::ActivateTimestep(): "
                 << msg.str() << std::endl;

          EXCEPTION2(NonCompliantFileException, "H5Part AcitvateTimestep",
                     msg.str());
        }
    } 
    else if (ts != activeTimeStep)
    {
        debug5 << "avtH5PartFileFormat::ActivateTimestep(): Activating time ";
        debug5 << "step " << ts << std::endl;

        if (H5PartSetStep(file, ts) != H5PART_SUCCESS)
        {
          std::ostringstream msg;
          msg << "Cannot activate time step " << ts << ".";

          debug1 << "avtH5PartFileFormat::ActivateTimestep(): "
                 << msg.str() << std::endl;

          EXCEPTION2(NonCompliantFileException, "H5Part AcitvateTimestep",
                     msg.str());
        }
    }

    // Update global time step information
    activeTimeStep = ts;

    visitTimer->StopTimer(t1, "H5PartFileFormat::ActivateTiumestep()");
}

// ****************************************************************************
//  Method: avtH5PartFileFormat::GetSubBlock
//
//  Purpose:
//      Compute extents of a sub bliock for a given processor id
//
//  Arguments:
//      gridDims      Number of nodes in x-/y-/z-direction
//
//  Output:
//      subBlockDims  Range of samples that this processor should read. Order
//                    i_min, i_max, j_min, j_max, k_min, k_max. Indices are
//                    inclusive, i.e., processor should also read samples i_min
//                    and i_max as well as all samples with indices in between.
//
//  Programmer: kurts and ghweber
//  Creation:   Tue Aug 28 17:35:50 PDT 2007
//
//  Modifications:
//    Gunther H. Weber, Thu Apr  1 15:43:05 PDT 2010
//    Use PAR_Rank() instead of block number as parameter. Fix problem
//    that cuased data to be only read partially if domain
//    decomposition was disabled.
//
// ****************************************************************************

void avtH5PartFileFormat::GetSubBlock(h5part_int64_t gridDims[3],
        h5part_int64_t subBlockDims[6])
{
    int t1 = visitTimer->StartTimer();
#ifdef PARALLEL
    if (enableDomainDecomposition)
    {
        int partitionAxis =0;
        int axis1 = 1;
        int axis2 = 2;

        if (gridDims[1] > gridDims[0] && gridDims[1] > gridDims[2])
        {
            partitionAxis = 1;
            axis1 = 0;
            axis2 = 2;
        }
        else if (gridDims[2] > gridDims[0] && gridDims[2] > gridDims[1])
        {
            axis1 = 0;
            axis2 = 1;
            partitionAxis = 2;
        }
        
        debug1 << "Grid dims are "
               << gridDims[0] << " " << gridDims[1] << " " << gridDims[2]
               << " partition axis is " << partitionAxis
               << " other axes are " << axis1 << " " << axis2 << std::endl;

        subBlockDims[2*axis1+0] = 0;
        subBlockDims[2*axis1+1] = gridDims[axis1] - 1;

        subBlockDims[2*axis2+0] = 0;
        subBlockDims[2*axis2+1] = gridDims[axis2] - 1;

        int numPartsForPartition = std::min(h5part_int64_t(PAR_Size()),
                                            gridDims[partitionAxis]);
        float slicesPerRank =
          float(gridDims[partitionAxis]) / float(numPartsForPartition);
        
        debug1 << "Slices per rank: " << slicesPerRank << std::endl;

        subBlockDims[2*partitionAxis+0] = int(PAR_Rank()*slicesPerRank);
        subBlockDims[2*partitionAxis+1] = int((PAR_Rank()+1)*slicesPerRank);
        
        if (PAR_Rank() == numPartsForPartition - 1)
        {
            // Make sure that the final slice is included and do not
            // include anything beyond.
            subBlockDims[2*partitionAxis+1] = gridDims[partitionAxis] - 1; 
        }
        
        debug1 << "Rank " << PAR_Rank() << " has block "
               << subBlockDims[0] << " " << subBlockDims[1] << " "
               << subBlockDims[2] << " " << subBlockDims[3] << " "
               << subBlockDims[4] << " " << subBlockDims[5] << std::endl;
    }
    else
    {
        debug1 << "Domain decomposition disabled. Returning entire block."
               << std::endl;
#endif
        subBlockDims[2*0+0] = 0;
        subBlockDims[2*0+1] = gridDims[0] - 1;

        subBlockDims[2*1+0] = 0;
        subBlockDims[2*1+1] = gridDims[1] - 1;

        subBlockDims[2*2+0] = 0;
        subBlockDims[2*2+1] = gridDims[2] - 1;

        debug1 << "avtH5PartFileFormat::GetSubBlock() returning block "
               << "comprising entire data set:"
               << subBlockDims[0] << " " << subBlockDims[1] << " "
               << subBlockDims[2] << " " << subBlockDims[3] << " "
               << subBlockDims[4] << " " << subBlockDims[5] << std::endl;
#ifdef PARALLEL
    }
#endif    

    visitTimer->StopTimer(t1, "H5PartFileFormat::GetSubBlock()");
}

// ****************************************************************************
// Method: avtH5PartFileFormat::GetDecomp
//
// Purpose: 
//   Distribute a set of objects evenly amongst all ranks
//
// Arguments:
//   nObjects   : The number of objects to be distributed
//   firstIndex : The index of the first object on this rank
//   lastIndex  : The index of the lastst object on this rank
//   total      : The total number of ojects on this rank
//
// Programmer: Allen Sanderson
// Creation:   Fri Jan 20 2017
//
// ****************************************************************************
void avtH5PartFileFormat::GetDecomp( h5part_int64_t nObjects,
                                     h5part_int64_t &firstIndex,
                                     h5part_int64_t &lastIndex,
                                     h5part_int64_t &total )
{
#ifdef PARALLEL
  h5part_int64_t rank   = PAR_Rank();
  h5part_int64_t nProcs = PAR_Size();

  h5part_int64_t nObjectsPerProc = (nObjects / nProcs);
  h5part_int64_t oneExtraUntil   = (nObjects % nProcs);
    
  if (rank < oneExtraUntil)
  {
    firstIndex = (rank  ) * (nObjectsPerProc+1);
    lastIndex  = (rank+1) * (nObjectsPerProc+1) - 1;

    total = lastIndex - firstIndex + 1;
  }
  else if( nObjectsPerProc ) 
  {
    firstIndex = (rank  ) * (nObjectsPerProc) + oneExtraUntil;
    lastIndex  = (rank+1) * (nObjectsPerProc) + oneExtraUntil - 1;

    total = lastIndex - firstIndex + 1;
  }
  else
  {
    firstIndex = -1;
    lastIndex  = -2;
    total = 0;
  }
#else   
  firstIndex = 0;
  lastIndex = nObjects-1;
  total = nObjects;
#endif
}

// ****************************************************************************
// Method: avtH5PartFileFormat::getVariablePathPrefix
//
// Purpose: 
//   Get the variablePathPrefix
//
// Programmer: Allen Sanderson
// Creation:   31 May 2017
//
// ****************************************************************************
const std::string
avtH5PartFileFormat::getVariablePathPrefix( int timestep )
{
  std::ostringstream prefix;

  prefix << variablePathPrefix << timestep;

  return prefix.str();
}
